#!/bin/env Rscript

library(mHMM)

argss <- commandArgs(trailingOnly = TRUE)

print(argss)

dat_file <- argss[2]
out_prefix <- argss[3]
out_name = paste(out_prefix, "cnv.mhmm.txt", sep=".")
c_const <- as.double(argss[4])
corr=T
if (argss[5]=='F') {
	corr =F
	out_name=paste(out_prefix, "no_corr.cnv.mhmm.txt", sep=".")
}

dat=read.table(dat_file, col.names=c('chr', 'pos', 'count1', 'count2'), header=T)

data_path=Sys.getenv('data_path')
chr_len=paste(data_path, 'b37_chr_ln_y.RData', sep='/')
load(chr_len)
# target_chr_len=b37_chr_ln[which(b37_chr_ln$chr %in% chr),]
target_chr_len=b37_chr_ln

dat$pos1=NA
for (i in 1:nrow(target_chr_len)) {
	dat$pos1[which(dat$chr==target_chr_len$chr[i])]=
	dat$pos[which(dat$chr==target_chr_len$chr[i])] + 
	sum(as.numeric(target_chr_len$len[0:(i-1)]))

}

sample_vec <- dat$count1
control_vec <- dat$count2
posi_vec <- dat$pos1

# sample_vec <- as.matrix(sample_vec, mode="numeric")
# sample_vec <- as.vector(sample_vec, mode="numeric")
# control_vec <- as.matrix(control_vec, mode="numeric")
# control_vec <- as.vector(control_vec, mode="numeric")
# posi_vec <- as.matrix(posi_vec, mode="numeric")
# posi_vec <- as.vector(posi_vec, mode="numeric")


# reference_reads_input <- control_vec
# sample_reads_input <- sample_vec
# loci_input <- posi_vec
kmeans1 <- TRUE

# change_point_refine <- TRUE

CNV_detection_result <- m_HMM_main(sample_vec, control_vec, posi_vec, kmeans1=kmeans1, C.const=c_const, change_point_refine = corr)

# filename <- paste(out_prefix0, "txt", sep=".")

# "CNV_detection_before_refine.RData"
# "CNV_detection_refine.RData"       
# "combined_result.RData"

if (argss[5]=='F') {
	write.table(CNV_detection_before_refine, out_name, quote = F, sep = "\t", row.names = F, col.names = F)
}else {
	write.table(CNV_detection_result, out_name, quote = F, sep = "\t", row.names = F, col.names = F)

}
source(paste(Sys.getenv('tools_path'), 'script/cnv_p_mhmm.r', sep='/'))
cnvpmhmm(out_name)



q()




